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Abstract 

A knowledge of the radial distributions of quarks inside hadrons could lead to 
a better understanding of the QCD description of these hadrons and possibly 
suggest forms for phenomenological models. As a step in this direction, in 
an earlier work, the charge (vector) and matter (scalar) radial distributions 
of heavy-light mesons were measured in the quenched approximation on a 
16 3 x 24 lattice with a lattice spacing of a ~ 0.17 fm, and a hopping parameter 
corresponding to a light quark mass about that of the strange quark. 

Here several improvements are now made: 1) The configurations are gen- 
erated using dynamical fermions with a ~ 0.14 fm; 2) Many more gauge 
configurations are included; 3) The distributions at many off-axis, in addition 
to on-axis, points are measured; 4) The data analysis is much more complete. 
In particular, distributions involving excited states are extracted. 

The exponential decay of the charge and matter distributions can be de- 
scribed by mesons of mass 0.9±0.1 and 1.5±0.1 GeV respectively — values 
that are consistent with those of vector and scalar gg-states calculated directly 
with the same lattice parameters. 
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I. INTRODUCTION 



In few- and many-body systems radial and momentum distributions often play an impor- 
tant role. For atomic and nuclear systems these distributions are, in many cases, calculated 
from a differential equation using an effective interparticle potential, where both the equa- 
tion and potential have some justification. However, for quark-gluon systems this approach 
is thought not to be applicable, even though the basic interaction — that of QCD — is 
exact and well known. Therefore, if — for example — transition rates between states in a 
heavy-light meson (Qq) are calculated, then the necessary radial wavefunctions are simply 
taken to have some convenient form, as in Ref . |TJ , or are calculated with a differential equa- 
tion and interquark interaction that are not well justified. It is an unusual situation, where 
one of the most fundamental systems, the hydrogen atom of quark physics, has interparticle 
correlations that are little understood. 

In an attempt to remedy partially this problem, the authors in Ref. 0] measured the 
charge and matter distributions in a heavy-light meson. More explicitly, the heavy-light 
meson was simplified to being an infinitely heavy quark (Q) and an antiquark (q) with a 
mass approximately equal to that of the strange quark. The physical meson nearest to 
this idealised meson is the _B S (5.37 GeV). Ref. was essentially a pilot calculation to test 
the feasibility of such distribution measurements and its successful outcome encouraged the 
authors to continue this line of research. In this paper, the same correlations are measured 
but with several major improvements compared with the earlier study: 

1) The gauge configurations are now calculated using two flavors of dynamical fermion 
and not in the quenched approximation as before. The actual parameters are those in 
Ref. ||, namely, (3 = 5.2, C$w — 1-76 for the clover action, a « 0.14 fm for the lattice 
spacing, k = 0.1395 for the hopping parameter and M PS /M V = 0.72 for the pseudoscalar 
meson/vector meson mass ratio, which corresponds to a quark mass somewhat heavier than 
the strange quark mass — Mps/My = 0.682 being the ratio corresponding to exactly the 
strange quark mass [|J. In comparison, the parameters for the quenched work of Ref. |2| 
were (3 = 5.7, C sw = 1-57, a « 0.17 fm, « = 0.14077 and M PS /M V = 0.65, corresponding 
to 10] a value of nig = 0.91 (2)m s . The spatial lattice size is 2.24 fm (compared to 2.72 fm in 
the previous quenched calculation). 

2) Many more gauge configurations are generated — 78 with dynamical fermions com- 
pared with the earlier 20 quenched configurations. 

3) Previously the densities were measured only at the 7 on-axis points r = 0, 1, . . . , 6. 
Now the densities at 20 points are measured — the same 7 on-axis with the addition of 
13 off-axis points. This permits a potentially much more detailed mapping of the density 
profiles. Furthermore, it opens up the possibility of testing the rotational invariance of 
these profiles by comparing them at r = 5 with that at (x = 3, y — 4) and also r = 3 
with (x = 2, y = 2, z — 1). This symbolic notation for an off-axis point will be used 
throughout the article. In practice, it includes the 8 possibilities (x = ±3, y = ±4) and 
(x = ±4, y = ±3), so that when combined with the directions along the 3 axes, this makes 
24 independent measurements for each symbolic (x, y), when x ^ y. In contrast, the on- 
axis cases have only 6 independent measurements for each r. This will be seen to improve 
significantly the statistics for the off-axis points compared with their neighboring on-axis 
points. This is in spite of the fact that an off-axis point requires a longer string of latticized 
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links, so that — being guided by strong coupling arguments — it should be more difficult 
to measure. 

4) Since the work in Ref. 0, the methods and our understanding of the data analysis 
have been developed. In particular, the interesting off-diagonal density terms are always 
allowed to vary and are no longer fixed to zero as was sometimes the case earlier. This 
now gives a better estimate of excited state effects. Also for the radial dependence of the 
density x a/3 (r), the use of a separable form y a (^)yp(r) is also employed. Here a, [3 are state 
indices with a = 1 being the ground state. This separable form is found to have some 
interesting features not found in the non-separable approach. Furthermore, if the densities 
were to be interpreted in terms of underlying wavefunctions ip a (r), then the separable form 
x a/3 (r) = ip a (r)ip/3(r) would be perfectly natural. However, it should be emphasised that 
such an interpretation can only have a phenomenological justification. 

In Section |H the two- and three-point correlation functions needed to extract the densities 
are briefly discussed — the reader being referred to Ref. |2j for more details. In Section [ft I 
the methods for analysing the basic lattice data are described. This results in values for the 
ground and some excited state energies and, in addition, radial distributions of the charge 
(vector) and matter (scalar) densities for these states are extracted. In Section |TV| these 
radial distributions are parametrized in terms of latticized forms of Yukawa, exponential 
and gaussian functions. In Section [V] a summary and some conclusions are made. 

II. THE CORRELATION FUNCTIONS C 2 AND C 3 (R) 

In this work the basic entities are the two- and three-point correlations [C2, C3], both 
of which are needed for measuring radial distributions. These are depicted in Figures [I] and 
H and are seen to be constructed from essentially two quantities — the heavy (static)-quark 
propagator Gq and the light-quark propagator G q . 

As discussed in detail in Ref. 0, when the heavy-quark propagates from site (x, t) to 
site (x', t + T), Gq can be expressed as 

G Q (x, t- x', t + T) = i(l + 74 )f/ Q (x,t,T)5 x , x ,, (I) 

where t/ Q (x, t, T) = Y[J=o Ui(x,t + i) is the gauge link product in the time direction. On 
the other hand, as the light-quark propagates from site i to site j, it can be schematically 
expressed as || 

G q = Gji = Qjl = ((QikfaYh) = (2) 

or as 

G' q = G' 3l = 75<(Q^k)#>75 = 75<^*>75- (3) 

Here the (pi are pseudo-fermions situated on the lattice sites i and ipi = Qik4>k, where Q is 
the Clover- Wilson-Dirac matrix, which is specified by C*sw an d the hopping parameter k. 

Knowing Gq and G q , the general form of a two-point correlation can be constructed from 
a heavy-quark propagating from site (x, t) to site (x', t + T) and a light-quark propagating 
from site (x', t + T) to site (x, t) as 
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C 2 (T) = Tr(rt G Q (x, t; x', t + T)T C ? (x', t + T; x, t)) 

= 2(Re [£/ Q [^*(x, t + T)0(x, t) + 0*(x, t + 7>(x, t)]l ). 



(4) 



Here T is the spin structure of the heavy-quark light-quark vertices at t and t + T. In this 
case T = 75, since we are only interested in pseudoscalar mesons such as the S-meson. For 
clarity, the Dirac indices have been omitted. 

Similarly, when the light-quark field is probed by an operator 0(r) at t — and the 
heavy-quark propagates from t = —t 2 to t — ti 



C 3 (-t 2 , t u r) = Tr(rt G Q (x, -t 2 ; x, t t )T G,(x, *i;x + r, 0) 6(r) G;(x + r, 0;x, -t 2 )). 



Here O = 74 for the charge distribution and 1 for the matter (scalar) density. 

The above has been written down for a single type of gauge field. However, the correla- 
tions can be greatly improved by fuzzing. In this case the basic links containing the gauge 
field have two fuzzings in addition to the original local field (L). In the standard notation 
of, for example Ref. ||, Fuzzl has 2 iterations and Fuzz2 a further 6 iterations i.e. 8 in all. 
These will be referred to as Fi and F 2 . In both cases, the factor multiplying the basic link 
is f p = 2.5 i.e. 

[Projected fuzzed link] = f p - [Straight link]+ [Sum of 4 spatial U-bends] 
with the quarks separated by a product of fuzzed links of length 1 lattice unit for Fuzzl and 
2 lattice units for Fuzz2, as discussed in Ref. @. With fuzzing included, C 2 and C3 are now 
3x3 matrices composed of matrix elements with the indices LL, LF 1; LF 2 , F]Fj, F!F 2 and 
F 2 F 2 . This means that S-wave excited state energies and properties can now be studied in 
addition to those of the ground state. 



There are several ways of analysing the above correlation functions C 2 and C3 in order 
to extract the quantities of interest i.e. energies and radial distributions. For a review of 
these methods see Ref. - with more details using the present notation being found in 
Ref. [Q. We now draw upon experiences learnt in that reference. 

Firstly, the two-point correlation data C 2 are analysed to give the energies (m a ) and 
eigenvectors (v) for the states of the Qg-system. These values of m a and v are then fixed 
when analysing the three-point correlation data C3 to give the charge and matter densities 



Consider the correlation function C 2 (T) as an n x n matrix — 3 x 3 in this case with 
the elements LL, LFi, . . ., F 2 F 2 . Each element C 2 ^j(T) is then fitted with the form 



(5) 



III. ANALYSIS 



x 




A. Analysis of the two-point correlation functions C 2 



C 2) y(T) w C 2 ^(T) = ]T i;?exp(-m a 7>; 



3 ' 



a 



(6) 
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where M 2 is the number of eigenvalues and mi is the ground state energy of the heavy-light 
meson. The values of m a and ■ are then determined by minimizing the difference between 
the C 2 data from the lattice and the form C 2 . The function actually minimized is the usual 



T 2 , 



x 2 = E E 

i>3 22,min 



C2,ij(T) — C2 ) ij(T) 
AC 2 , i3 -(T) 



(7) 



where AC 2) y(T) is the statistical error on C 2)i j(T) and T 2imin , T 2max are the minimum and 
maximum values of T 2 used in the fit. Here T 2 max is fixed at 11 — the maximum possible 
value when applying the Maximum Variance Reduction method on a 16 3 x 24 lattice — see 
Ref. ||. The minimization is carried out using the program package Minuit — the Migrad 
option being the most successful, since this enables errors on the varied parameters m a and 
vfj to be determined. It was unnecessary to use the Minos option, since the errors were 
found to be quite symmetrical. 

Before the above minimization is started, two parameters need to be decided upon — M 2 
and T 2im i n . The criteria for this choice are that the x 2 / n dof must be small — to be satisfied 
that the fit to the data is good — and that T 2jmin is appropriate for future use of the m a and 
vfj in analysing the C 3 data. The latter data covers a range from T 3 = 4 to 10, where each 
C 3 is constructed from essentially the product of two C 2 's i.e. C 3 (T 3 = 10) requires a good 
fit to C 2 (T 2 = 5), whereas C 3 (T 3 = 4) would require a good fit to C 2 (T 2 = 2). However, 
it must be kept in mind that ultimately, when extracting densities from C 3 (T 3 ), it is the 
larger values of T 3 — where convergence hopefully occurs — that are of more importance. 
In practice, this means T 3 min « 8 is a reasonable compromise, so that the C 2 {T 2 ) fit should 
concentrate on values of T 2 > 3. 

In Ref. [[|, M 2 = 2,3 and T 2 min = 3,4 were considered. However, the M 2 = 2 and 
^2,min = 3 fits resulted in considerably larger xV n dof and were not used in the subsequent 
C 3 analysis. Here the situation is similar with now M 2 = 3,4 and T 2im i n = 3,4 giving the 
best and most appropriate fits to C 2 . The outcome is shown in Table |, where we present 
4 cases. In most of this paper we will concentrate on Case B, since this has both a good 
X 2 '/ n dof of 0.16 and sufficiently small errors on the state energies am a . In contrast, Case A 
with M 2 = 3 has a large x 2 / n dof and Case C large errors on the am a . 

The table also shows — with Case Q — the earlier best fit in Ref. || to the 20 quenched 
configurations with (3 = 5.7. When comparing the am a from the four cases, two points must 
be kept in mind: 

1) Only differences of the am 's have a meaning, since the lattice simulation generates 
different self-energies to the quarks in Case Q versus Cases A, B, C. 

2) The table shows am a , where the lattice spacing a is ~ 0.17 fm for Case Q and ~ 0.14 fm 
for Cases A, B, C. 

Removing these two effects results in the Am^ = amp — am a at the bottom of Table |. 
There it is seen that, within the error bars, both Am 2 i and Am 31 are unchanged in going 
from the 20 quenched configurations — after being scaled by the ratio of the lattice spacings 
0.14/0.17 — to the preferred unquenched Case B. The best estimates are Am 2 i = 0.33(1), 
Am 3 i = 0.80(4) and Am^i = 1.05(6). This value of Am 2 i is also the same as that obtained 
in Ref. [H], when the latter is also scaled by the lattice spacing ratio 0.14/0.17. However, it 
should be pointed out that in Ref. the quenched approximation was used with the same 
parameters as in Ref. || but on a 12 3 spatial lattice. 
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In Table |, for Case B, there are 54 bits of data fitted with 16 parameters. However, 
it should be added that the Parameter Correlation Coefficient from the Migrad algorithm 
of Minuit indicates that several of these parameters are correlated. For example, there is a 
strong correlation between ami and and, if i>£ is fixed as = 0.45046 • ami, then Migrad 
finds the same solution. In this way 9 of the 16 parameters can be effectively removed to 
leave only 7 — the 4 rrii and v F F , vp 2 - - but with the same final solution. Therefore, 
the n2,dof increases from 38 to 47 and the x 2 / n 2,dof decreases to 0.13. Such small values of 
X 2 / n 2,dof are not surprising, since also correlations within the data are expected. These are 
of three forms: 

1) For a given Cy(T), there may be correlations between the 78 gauge configurations. How- 
ever, this has been greatly reduced selecting the 78 configurations used to be separated by 40 
trajectories. This feature can be checked by calculating the autocorrelation function (/) be- 
tween the N=78 configurations to give an effective number of configurations A^gective — NJ f. 
Here it is found that / ~ 1 — indicating a negligible correlation. 

2) For a given gauge configuration, there are correlations between the different T-values for 
both C*2(T) and C 3 (T). This has been studied in, for example, Ref. ||, where the general 
conclusion was that consideration of these types of correlation resulted in essentially the 
same final eigenvalues but, in some cases, with somewhat increased error bars. 

3) There could be correlations between the different values of r. Hopefully, our employment 
of various analysis schemes — such as the use of several different values of T min and also 
the non-separable versus separable forms discussed below — will minimize the uncertainty 
in this particular correlation. 

B. Analysis of the three-point correlation function for radial distributions 

The analysis of the three-point correlation functions C${®,T, r) is performed using a 
generalisation of the one for C2 in Eq. |6], namely, 

M 3 M s 

C 3)ii (e,T,r) w C 3iij (e,T,r) = £ £ ^expf-m^a^r) exp[-m p (T - t^vf. (8) 

a=l (3=1 

The m a and v-vectors are those obtained by minimizing the C2 in Eq. |6] and, for each 
value of r, the x al3 (r) are varied to ensure a good fit to Cs t ij(Q,T,r) by the model form 
C3,ij{Q,T,r). As for C2 a decision must be made on the maximum and minimum values 
2"3,max) ^3,min of T 3 . Now T 3 max is fixed at 10, since the signal to noise ratio was too large 
for larger values of T. As discussed above, when deciding on T 2imin , we consider T 3jmin « 6-9 
to be a reasonable range. 

Two forms of x a "{r) are used here: 

1) A non-separable (NS) form, where each x a/3 (r) is treated as a single entity. Here we 
take M 3 = M 2 = 4. However, for Migrad to converge to a reasonable solution, of the 10 
possible values of x a "(r) for a given value of r, only 7 are varied — the 4 x aa (r) and the 
3 x la (r) with a / 1. The other x al3 (r) are fixed to be zero. As discussed above for C2, 
Migrad indicates correlations between these 7 variables. However, taking into account these 
correlations results in the same solution. 

2) A separable (S) form x a ^(r) = y a {r)y/3(r). Here, we take M 3 = 3 to give only three free 
parameters for each value of r — yi{r), y2{r) and y 3 (r). 
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It is seen from Section |I| that this second form appears to be a more natural parametrization, 
since schematically we can write C 3 as {^$Gq(— t 2 — > ti)j 5 G q (ti — > 0)QG q (0 — > —t 2 )) in a 
form containing two light-quark propagators. Because these propagators act over different 
time intervals, G 3 ^ 3 - might — to some extent — be separable in terms of the indices a and 
(3. As shown in the Appendix, when the separable form x a/3 (r) = y a {r-)yp{r) is used in Eq. |8| 
the ratio 



Rij(T) = In 



C 3)i j (T 



(9) 



has, when T is odd, the symmetry 

R l3 {T) = R l3 {T + l). (10) 

The surprising point is that this same symmetry is also found, to a greater or lesser extent, 
in the basic data. This is shown in Table |H] f° r the dominant correlation C 3 ^ 1 ^ 1 at the small 
interquark distances of r = 1 and the off-axis points (x = 1, y = 1) and (x — 1, y — 1, z — 1) 
the symmetry is very clear. However, as r increases the symmetry essentially disappears - 
as seen in the last column for (x = 3, y = 4). 

In Table [HJ results are given for both the NS and S forms of x a @(r) and for different 



choices of T 2imin , T 3imin — with the most representative solution being S(3, 8), the separable 
form with T 2>min = 3 and T 3imin = 8. The other choices give support to this solution and 
indicate the possible systematic error. In the penultimate column is given, in our opinion, 
the best overall estimate of the ground state charge density with error. The last column 
shows estimates using Case A in Table |. As shown in Ref. |J, in this case the matrix 
is square and so can be inverted to give the matrix w". Estimates C^ a p (T) can then be 
written down directly as 

C 3 , a p(T) = u?C 3tij (T) U P. (11) 

In this case 

1 (r> = j"" w^m, ■ (12) 

Unfortunately, the extraction of the asymptotic T — > oo limit is somewhat subjective and 
gives the estimates in the last column in Table [TTT| . In all the cases listed, within errors 
these agree with the previous column. However, this approach did show that the data from 
some of the larger values of (x, y) were not good and so these are dropped in the subsequent 
discussion. Also, for reasons to be discussed in Section [IV E| , neither the (2,2,1) data nor 



the weighted average of (3,0,0)/(2,2,l) are included in the following analysis. 



In Table [TV] similar estimates are given for charge densities involving excited states. Also 
this table contains a summary of the x ll (r) and x 12 (r) matter radial distributions, which 
were extracted using both the NS and S forms for different choices of T 2imin , T 3imin — just 
as in the charge case. However, these signals are somewhat weaker than for the charge, so 
that no meaningful matter distributions could be extracted for r > 4a. 

In Figure ^| the best estimates of the charge and matter distributions from Tables |Hj 



and IV are compared. To guide the eye we also show lines depicting the lattice exponential 
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fits to be discussed later in Section [IV A| . Here it is clearly seen that the range of the 
charge distribution is longer than that of the matter distribution. Furthermore, this figure 
also contains the charge and matter densities obtained with the quenched approximation in 
Ref. @. For this comparison, the results of Ref. are scaled from lattice spacing a o u to 
a o.n by simply r .i7 — > pr .i7 and x u (r ,u) — > p~ 3 x 11 (r n), where p = 0.17/0.14. Here it is 
seen that the present results using dynamical fermions are indistinguishable from those using 
the quenched approximation. It will be seen later in the last two columns of Table IVJ, that 
this near equality is also reflected in the charge sum rules with both giving x 11 = 1.4(1). 
As discussed in Refs. || and |TD|, this is of particular interest in the matter case, since 



there disconnected contributions could arise that are dependent on the quenched versus 
unquenched. Any difference would then be due to the effect of the quark condensate. Clearly, 
with the present data no such effect can be detected. However, it must be remembered that 
here the sea quarks have the same mass as the valence quarks i.e. about that of the strange 
quark. It is possible that using sea quarks with u, d masses the above conclusion would be 
different. This observation that full QCD and the quenched approximation do not differ 
significantly was has been seen many times before. 

In Figure ^a), for the charge density ratio x 12 /x 11 , only the errors for the separable 
analysis with T 2)m i n = 3 and T 3jmin = 8 are shown, since the other analyses have similar 
errors. There a distinct node is seen at about 2.2 lattice spacings i.e. at ~ 0.3 fm. Such a 
node is natural for x 12 , since it involves the excited S-wave state. On the other hand, for 
the matter density — as seen in Figure ^c) — the node is near to 1.5a ~ 0.2 fm. Figure |]b) 
shows the various analyses for the charge density ratio x 13 /x u . Here the node structure is 
less clear. One node is seen at about 2.8 lattice spacings i.e. at ~ 0.4 fm. But a second 
possible node at about 0.6 lattice spacings, i.e. at ps 0.1 fm, depends on the one value of 
x 13 at r = 0. However, we have no reason to suspect that this is purely a lattice artefact. 
Furthermore, for second excited S-wave states a second node is not unexpected. Similar 
comments hold for the matter density ratio x 13 /x u in Figure |]d). The node structure of 
x 22 — the charge density of the first excited state — is not at all clear. If x 22 is expressed 
in the separable form 2/2( r )2/2(' r ) ? then the zero that should appear at about 2a is not seen 
very distinctly in comparison with that seen in Figures ||a) or c) for x l2 /x u . 

The above figures show directly the various charge densities x a ^(r). However, it is 
also of interest to see the structure of the individual terms y a {i r ) in the separable form 
r) = ya{f)yp{f)- These are shown in Figure |5] for both the charge and the matter. 



x 



a/3. 



Figure |5|a) shows clearly that y\(r) for the charge has a significantly longer range than for 
the matter. Also as seen in Figure |5|b) both of the t^O") exhibit a distinct node and are 
responsible for the nodes in the separable form of the density x 12 already seen in Figure Hk). 
We do not plot y 3 (r), since the signal/error ratio is too small. 

As discussed at the end of the Introduction, the y a (r) can possibly be interpreted as 
wave functions for the state a. However, there are other radial distributions associated with 
the Qq system that can also be interpreted as wave functions. These are the Bethe-Salpeter 
wavefunctions [if a (r)] discussed in Ref. |J. They were extracted by assuming the hadronic 
operators C^aatji, ?"2, T) to be of the form w Q ,(ri)w a (r 2 ) exp(— m a T), where the sink and 
source operators are of spatial size r% and r 2 . In Figure [| a comparison is made between the 
above values of yi(r) and j/ 2 (r) and the corresponding results form Ref. || for u>i(r) and 
w 2 (r), where the latter have been normalised so that u>i(0) = yi(0) and w 2 (0) = 2/2(0) and 
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the values of r scaled. Even though they do bear some similarities, it should be added that 
there are several reasons why these two types of wave function should not agree in detail with 
each other. In particular, the [u> Q (r)] 2 cannot be identified as a charge or matter distribution. 
In addition, they were found by using an explicit fuzzed path between Q and q and so are 
dependent on the fuzzing prescription, whereas the y a are defined in a path-independent 
way. Furthermore, in Ref. || the w a (r) were extracted in the quenched approximation on 
a 12 3 x 24 lattice with the parameters j3 = 5.7, C*sw = 1-57 and k = 0.13843 — the latter 
corresponding to a light-quark mass of = 1.77(4)m s i.e. about two strange quark masses. 



C. Charge sum rule 

In addition to measuring C^{r) for various values of r, the correlation where r is summed 
over the whole lattice is also obtained. This leads to the charge sum rule as discussed in 
Ref. M. The actual values of this sum rule are extracted using Eq. ^ where the x a/3 are 
now independent of r. The outcome — as shown in Table [V] - is that x 11 is ~ 1-3(1) 
consistent with the earlier quenched result. The fact that x 11 is not unity — as expected 
in the continuum limit — can be qualitatively understood by introducing a renormalisation 
factor of ~ 1/1.3 ~ 0.8 into the basic 74 vertex used to measure the charge density. Such a 



factor of this magnitude is reasonable as shown in Ref. [II . 

It is also reassuring that the x a @ with a 7^ (3 are, in general, consistent with zero - 
as expected in the continuum limit. However, the interpretation of x 22 is less clear. The 
non-separable (NS) case gives x 22 ~ 1-0(4) — again a reasonable value in the continuum 
limit — whereas the separable case (S) yields x 22 ~ 0.0. This suggests that the separable 
approximation may be less appropriate for excited states. 

In Table [V] we also show the matter sum rule. These have a somewhat wider spread 
of values with 0.9(1) being a reasonable compromise — a number that is about twice the 
estimate of 0.38(15) for the quenched calculation of Ref. f2|. Perhaps this is an indication 
- unlike the matter radial distributions in Figure [3] — that the quenched and unquenched 
results can differ even with the present sea quark masses of about the strange quark. How- 
ever, we do not have the data to cross check with Refs. || and | ID| , which advocate the 
existence of such a difference for the matter sum rule. 



IV. FORM OF RADIAL DEPENDENCE 



The results in Table [TTJ are presented as simply a series of numbers for each value of r 
or (x,y,z). However, it would be more convenient and perhaps illuminating, if they could 
be parametrized in some simple way. This can be done either in coordinate or momentum 
space — the topics of the next subsections. 



A. Coordinate space fit 



Here we assume that the radial dependences of the lattice data in Table [ITT] can be 
represented in terms of exponential (E), Yukawa (Y) or gaussian (G) functions. Then the 
strategy is to first fit the data at the largest values of r with a single form in order to 
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parametrize the longest ranged part of the data, where it is expected that lattice artefacts 
will be less. This range is then used as a starting point, when the data at all values of r are 
fitted by adding a second form. 

The reason for using exponential and Yukawa radial functions is that they arise naturally 
as propagators in quantum field theory — usually in their momentum space form (q 2 +m 2 ) . 
However, if we go away from quantum field theory and attempt to understand the radial 
dependences in terms of wavefunctions from, for example, the Dirac equation, then gaussian 
forms can then arise. 



1. Exponential, Yukawa and gaussian fits to the charge and matter densities. 

Assuming simply an exponential a E exp[— r/r E ], Yukawa a Y exp[— r/r Y ]/r or gaussian 
a G exp[— (r/r G ) 2 ] form and, if all 18 pieces of data are used (excluding that for r = in 
the Yukawa case), then the x 2 / n dof are greater than unity — being 1.7, 2.3 and 5.1 for E, 
Y, G respectively. However, these x 2 / n do{ can be reduced significantly by adding a second 
exponential, Yukawa or gaussian form with the first form having its amplitude a F and range 
r F guided by those values needed to fit the data with r > 3. When these combinations are 
used to fit all the data (excluding that for r = in the Yukawa case) we get the results in 
the columns marked as 2E, 2Y and 2G of Table [V|. We see that the three radial integrals 
J F are comparable at w 1-4(1) and in reasonable agreement with those obtained earlier by 
directly summing over the whole lattice and shown in Table 0. These will be discussed later 



in Section [IV D| , where a detailed analysis of the charge sum rule will be made. 

The various merits and defects of the E, Y and G forms will be discussed in more detail 
later in the context of the separable approximation for the density. At this stage all three 
forms are equally reasonable giving x 2 / n dof that are ~ 1.4. 

In Table |V11| are shown the corresponding fits to the matter density from column 5 of 
Table |V|. Here it is seen that for 2E and 2Y the x 2 / n do{ are ~ 1.7 — larger than in the 
charge case. In comparison, for 2G the situation is reversed with a x 2 / n dof = 1.11. 



2. Lattice exponential, Yukawa and gaussian fits to the charge and matter densities. 



As seen in Fig. |3|, the data is far from being a smooth function of r — indicating lattice 
artefact effects. Therefore, in addition to fitting the data with pure exponential, Yukawa 
and gaussian forms, lattice versions of these (LE, LY, LG) are also used. In Ref. |12| the 
lattice form of the Coulomb function (1/r) is written as 



7r ^ cos(r.q) 



LC 



aL 3 



E 



D 



(13) 



Here L is the lattice size along one axis and D = Z)f =1 sin 2 (agj/2), where aqi = 
0, ^r, --- 27r (^~ 1 ) ; q ^ o. In this subsection, for clarity, the lattice spacing a is shown 
explicitly. For the above Yukawa form, Eq. O is easily generalised to 



exp(— r/r ) 



7T 



cos(r.q) 



J LY 



aL 3 ^ D + 0.25[a/r LY } 2 ' 



(14) 
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However, now the point q = can be included in the sum, since it is no longer a singularity 
- provided l/r LY 7^ 0. Similarly, lattice forms of the exponential and gaussian can be 
obtained by simply replacing in the usual Fourier transform the q 2 factors by their lattice 
equivalent \ Yh=i sin 2 (agj/2) = ^ . This results in 



exp(— r/r 



E 



cos(r.q) 



le 2r LE L 3 Y \ D + 0.25(a/r 



LE^212 ' 



(15) 



exp[-(r/r LG ) 21 



LG 



r LG y/TT 

aL 



cos(r.q) exp[-(r LG /a) 2 D}. 



(16) 



When using these lattice forms, columns 2LY, 2LE and 2LG in Tables |V| and [VI I| for the 
charge and matter respectively show the results corresponding to the usual forms 2Y, 2E 
and 2G. There it is seen that the values of r$ for all three lattice forms 2LE, 2LY and 2LG 
are quite similar to their 2E, 2Y and 2G counterparts. The outcome is that the lattice 
forms are able to reproduce some of the structure especially near r = 2 and 3 — as seen 
in Figures || and |6|. This is also reflected in the xV n dof's being now less than unity for the 
charge case and < 1.2 for the matter. At first sight the two parametrizations 2LE and 2LY 



look very different, since Eqs. |1J and [15] have such dissimilar forms. In addition, in 2LE 
the two terms add up, whereas in 2LY they cancel, since Oq LY and a 2LY have opposite signs 
in order to dampen the 1/r effect at small r. However, later we shall see that, in practice, 
2LE and 2LY behave in very similar ways — with little numerical preference for one over 
the other. 

The conclusion from this subsection is that all three parametrizations 2LE, 2LY and 
2LG are acceptable, since each can fit all the lattice data with x 2 / n dof's that comparable to 
unity. 



3. Exponential, Yukawa and gaussian fits to the separable form of the charge and matter densities. 

The above has concentrated on fitting directly the ground state charge density x n (r) 
in Eq. ||. However, in section |ll IB] a second and possibly more natural parametrization - 
a separable form x a/3 (r) = y a (r)yp(r) — was introduced. This resulted in the data shown 
in columns S of Table |T|. Here, we consider that the j/i(r) are simply y / x 11 (r), where the 
latter are the Best Estimate values in Table TT|. The corresponding fits 2ES, 2YS and 2GS, 
with the above forms are shown in Tables |VTII| and [IX|. For reasons to be seen later, we only 
show the matter results for 2YS. Several points can be seen: 

1) As expected, the ES and GS fits are now simply related to those for the density — namely 
the (Zq ~ (<2q S ) 2 , Tq S « 2r E and r GS « V2r G . Therefore, since little is added by the E and 
G fits, we do not show them for the matter case. 

2) When fitting all the data, the 2ES, 2YS and 2GS fits are comparable but all with xV^dof's 
that are greater than unity. When data points at small r are systematically removed, then 
the x 2 / n dof eventually goes under unity. In the charge case for 2ES, 2YS and 2GS this 
happens at r > 2.24, 2.24 and 2.0 respectively. However, the fact that the three x 2 / n dot 
are large is partially due to the lattice artefacts already seen when dealing directly with the 
density in the previous subsection. 
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4- Lattice exponential, Yukawa and gaussian fits to the separable form of the charge and matter 

densities. 



Here, the yi(r) are fitted using the lattice exponential, Yukawa and gaussian forms in 
Eqs. [14], [15] and [16] and the outcome is shown in Tables |VIII| and [D|. The fits 2LES, 2LYS, 
2LGS are slightly worse than before — but still acceptable. In 2LYS the parameters ao and 
ai are highly correlated with a\ ~ — ao — again to dampen the 1/r effect for small r. In 
the limit for small r the 1/r would cancel exactly to leave a function more like an 

exponential. 



The conclusion to be drawn from this Section is that the lattice data extracted in Sec- 
tion III can be well fitted by any of the three lattice forms in Eqs. 14, 15 or |16| — numerically 
none of them is superior and also none can be rejected. This statement applies not only 



to the direct parametrizations of the density as in Section [IV A 2| but also to the separable 
form in Section II V A 4. 



B. Momentum space fit 

Often it is more convenient to view data in momentum space by making the transforma- 



tion 



x af} (k) = VJcos(k.r)x a/3 (r), 



(17) 



where the r summation should be over the whole 3-dimensional L 3 lattice. This would mean, 
in the present calculation, summing x, y and z over the ranges —7 to +8. If the lattice 
Yukawa, exponential and gaussian expressions in Eqs. [14], [15| and |l^, written as 

5>os(r.q)F(q), (18) 



x r 



7T 



are now used to parametrize x(r), then we simply get 

x(k) = ?rF(k). 



(19) 



Here use has been made of the identity YJ r cos[r.(k — q)] = L 3 0k, q . From Eq. [17] we see that 
the sum rule discussed earlier is now directly x(k = 0) = 7rF(0) = I F , where the I F are 
defined in the caption of Table |VI. The other fourier components with k / are simply 
vrF(k). 



1. The inclusion of lattice data directly into Eq. 11 



The result in Eq. [19] is not surprising, since it is simply the Fourier transform of a Fourier 
transform. However, if — where available — the lattice data of Table [TTT] are used for the 
x a P(r) in Eq. 17 and the form in Eq. [18] only used for the missing densities, then the result 
should be an improved estimate of the x a/3 (k). A measure of this can be obtained by fitting 
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this improved estimate with the forms F(k, a^r^) in Eqs. |14|, [15] and [16|, but where the 
parameters tuned using Minuit and are not the ones appearing in Tables [VI] and 



As we shall see in the next subsection, the fits of main interest are 2LY, 2LE and 2LYS, 
since they possibly have a physical interpretation. However, this fitting procedure presents 
some problems, since we are attempting, with only four parameters (ai=o,i, rj = o,i), to fit 165 
numbers — the x a/3 (k) with hi — 0, . . . , j n particular, if — as a stability check - 

an analytic form of F(k, aj,r«) is analysed using the same form, then the Migrad algorithm 
of Minuit does indeed give exactly the correct values of the ai,rt but with large error bars. 
Another problem is the choice of function to be minimized. We consider two options: 

(i) iVi = [x FT (k) - x mode i(k)] 2 and 

(ii) N 2 = [x FT (k)/x model (k) - l] 2 , 

where the x F T(k) are the Fourier transforms (with or without direct data inclusion) defined 
in Eq. [17| and x mo dci(k) are the 7rF(k) in Eq. [l^. Since the x F T(k) decrease rapidly as k 
increases, Option (i) emphasizes the smaller values of k and is appropriate for extracting ro, 
which is the longer range. In principle, Option (ii) is better since it should give an overall fit 
to the x F x(k). But in practice, it tends to be unstable yielding either unacceptable solutions 
or very large error bars. 

If the 2LY fit to the charge density in Table [V]] is analysed with a Yukawa form, then 
Options (i) and (ii) give poles at r = 1.78(24) and 1.78(9) respectively — see Table This 
is to be compared with the input value of 1.78(7) in Table [VI]. A similar strategy can be 
applied to the improved x u (k) generated from Eq. [T7|, where x u (r) now contains directly 
lattice data wherever possible. In this case, Options (i) and (ii) yield ro = 1.78(20) and 
3.3(1.0) i.e. the result from Option (i) is not distinguishable from using the fitted expression 
in Eq. |14] for all values of r, whereas that from Option (ii) is unstable. Similarly, if the 2LE 
fit in Table [VT] is analysed with an exponential form, then Options (i) and (ii) give poles 
at r = 1.38(22) and 1.38(15) respectively. This is to be compared with the input value of 
1.37(6) in Table [Vl|. For the improved form of the density the two Options yield 1.39(23) 
and 1.41(16). Again Option (i) is indistinguishable from using only the fit values of a; 11 (r). 
However, now option (ii) is more stable than before. 

In addition to analysing Yukawa forms with Yukawas or exponential forms with expo- 
nentials other combinations are possible. These are summarised, alongwith the above, in 
Table |X[ Here the main interest is that, when analysing the exponential forms with a 
Yukawa, the same value of r = 1.75(24) emerges. 

The conclusion to be drawn from Table [X| is that the parametrizations in Tables |V1] and 
Vlll| are so good that any improvements on r due to the inclusion of explicit lattice data 



cannot be detected. Therefore, in the following discussions the Fourier transforms based 
purely on Tables [V| and Vlll will, in general, be used. 



C. Possible interpretations of the above fits. 

In the above, the use of the various forms E, Y, . . . , LYS, LGS is considered as a 
purely numerical exercise. However, one can also ask about any theoretical interpretation 
of, or preference for, one form over the others. As we shall see below, in some cases, this 
is best discussed in momentum space. Also, since the inclusion of the original data has 
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little effect, the fourier transform of the charge density should be well described by simply 
7rF(k). Therefore, the discussion below focusses on the interpretation of the different forms 
of 7rF(k). Here several possibilities are suggested: 

1) Y and LY. These forms can be directly identified as the propagators of a single 
particle. In the Y form — aj exp[— r/rj]/r — the masses of the propagating particles are 
simply 1/rJ . Also for LY by writing the denominator of Eqs. [14] in the continuum limit as 
2j- + 0.25[a/rf Y } 2 , we see that these masses can also be identified as l/rf Y . For the charge 
case these will be vector particles, whereas for the matter case they are scalars. From the fits 
to the charge distributions in Table VI , the vector masses that emerge are — in lattice units 
of about 1.4 GeV — am^ 2LY = 

0.56 ±0.02 and 

am i2LY = 0.89±0.11 with the corresponding 
scalar masses from the matter distribution in Table |VI]] being am^ 2LY = 1-09 ^ 0-06 and 
am\ 2ljY = 1.4 ± 0.2. 

2) E and LE. Looking at Fig. ^| for the three-point correlations, a cut in the T-direction 
intersects the two light-quark propagators G q (t 1 — > 0) and G q (0 — > —t 2 ). In comparison, 
the lattice exponential form in Eq. [1^ also contains a product of two propagator-like terms 
1/[D + 0.25(a/r E ) 2 ]. This then suggests that the exponential form may be interpreted as 
the product of two non-interacting quark propagators. Assuming that the momentum of 
the probe is divided equally between the two propagators, then the appropriate momentum 
transfer in each propagator becomes q/2. Now, when going to the continuum limit, the 
denominators in Eqs. |15| give the masses of the propagating particles as l/2rf E i.e. in 
the charge case am 0j2LE = 0.36 ± 0.02 and am 12 LE ~ 4.5. Therefore, one interpretation 
of ^o,2LE is that this is the mass of a constituent quark used in the naive quark model 
description of the meson as simply two non-interacting quarks. This would give a vector 
mass of amQ 2 le = 0-72 ± 0.04. A similar interpretation can be made for the matter case to 
give a scalar meson of mass am^ 2LE = 1.07 ± 0.05. 

It should be added that in our earlier work in Ref. |2[ using the quenched approximation, 
the data were so sparse that an overall fit with only a single exponential was attempted 
over a limited range of r values — the overall fit with a single Yukawa being much worse. 
Therefore, to compare with the above values of am v,s , a single Yukawa fit to the data at 
the largest values of r for which the data were still reliable was carried out. For case 3 in 
Ref. 0, the charge density data at r = 3 and 4 gave am v =0.6(1) — scaled to a = 0.14 fm. 
Similarly, the matter density data at r = 2 and 3 gave am s =1.0(2). These values are not 
significantly different to the present estimates. 

In the above, the masses have been extracted by a somewhat tortuous argument. How- 
ever, in the literature there have been direct calculations of the energies of these qq states 
using the same lattice parameters and lattice size as those employed here. In Ref. || they 
got amp = 0.785 ± 0.009 and in Ref. [0] am^ = 1.18 ± 0.08. These numbers are consistent 
with our above estimates from the 2LE fit i.e. 0.72 ±0.04 and 1.07 ±0.05 respectively. Also 
from the 2LY fit, the scalar mass of amQ 2hY = 1-09 ^ 0-0^ * s consistent with this value. 
But the vector mass of cm^LY = 0.56 ± 0.02 is somewhat too small. This difference can 
be considered as a measure of the systematic error on m[j and suggests that the present 
data do not extend to sufficiently large values of r for a reliable estimate to be made of the 
asymptotic form. 

Unfortunately, it is not straightforward to identify the above particles and their masses 
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directly with physical particles, since we use light quarks that have the isotopic spin proper- 
ties of u, (i-quarks but with masses about that of the strange quark. In addition, we do not 
calculate mass or density contributions that arise from disconnected correlations. For vector 



mesons, the latter has been shown in Ref. [13[ to be only a small effect i.e. the OZI rule is 



justified in this case. On the other hand, this appears not to be so for scalar mesons. Evenso, 
in the vector case it is not reasonable to identify the above range of values itiq = 0.9 ±0.1 
and m\ = 1.4 ± 0.3 GeV directly with the isovector p(0.77) and the radial excited p(1.45) 



from Ref |14j], since in the quark model, the latter are constructed from u, d quarks with the 
correct mass — a value much less than the strange quark mass used here. However, since 
the OZI rule is a very good approximation for vector mesons, the additional mass of a state 
with strange quarks can be taken from the 0(1.02) meson — a value somewhat larger than 
our estimate of m^J = 0.9 ±0.1 GeV. For the scalar mesons, the comparison with experiment 



is even more indirect. When comparing with our results, from Ref. [|14[] the appropriate 



states would be the a (0.98) and a (1.45), since our neglect of disconnected correlations 



effectively results in an isovector operator. However, as recently discussed in Ref. [15| , these 
states probably have a complicated structure being mainly (qq)i(qq)3, in S-waves at short 
distances, with some qq in P-waves. But further out these rearrange into {qq)i{qq)i and 
finally emerge as meson-meson states. 

In fact this complication arises also in lattice calculations, since the present lattice pa- 
rameters predict a pseudo-scalar meson (the "pion") with mass am ps =0.564 — see Table 
8 in Ref. || . This means that in Fig. a cut in the T-direction can intersect four light- 
quark propagators — a state that can be interpreted as the exchange of two pseudo-scalar 
mesons. In the above scalar meson the qq are in a relative P-wave, so that it can couple to 
a two-meson state, where the mesons are in a relative S-wave. On the lattice such a two 
pseudo-scalar meson state would simply have a mass of 2am ps =1.13. This means that the 
scalar meson — calculated directly as a qq state with mass am s Q = 1.18 ± 0.08 in Ref. |7j 
- is essentially degenerate with the two pseudo-scalar meson state and suggests that our 
estimate of, say, amg 2 ly = 1-09 ^ 0-07 a l so contains this two-meson effect. It should be 
added that this problem does not arise for the vector meson, since there the qq are in a 
relative S-wave, so that the corresponding two-meson state has a relative P-wave leading to 
an energy considerably larger than 2am ps =1.13. This means that the structure of the vector 
meson generated here is expected to be mainly qq with little mixing with the two-meson 
state. 

3) YS and LYS. The above considers directly the charge density. However, a similar 
analysis can be carried out on the yi(r), when the charge density is expressed in the separable 
form x n (r) = y\{r)y\{r) . Since the single form fits ES, GS, LES and LGS are given directly 
by their non-separable counterparts, as is clear by the relationships a F ~ (a FS ) 2 , r FS « 2r F 
and r$ s rs V^r^, they would add nothing new. The same is essentially true of the 2ES, 
2GS, 2LES and 2LGS fits in Table Will, since the values of the 7- FSy/LFS and r 2FS//2LFS are 



very similar. However, if the y%(r) are thought of as one of the two propagators in Figure 



then it is now appropriate to identify yi(r) with the YS and LYS forms. From Table [VTlI 
we get that m 0i2L YS = 0.49 ± 0.02 and m 1>2LY s = 0.58 ± 0.02 GeV. 

This value of mo,2LYS is consistent with the earlier mo,2LE = 0.51 ± 0.03 GeV and sup- 
ports the identification that these two very different ways of analysing the data are indeed 
extracting the propagator of the same "particle" and, as said above, it is tempting to identify 
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this "particle" with a constituent quark. Since the density in the 2LE case can be expressed 
schematically in terms of two Yukawas (Yj) as YqYq + YiYi compared with the corresponding 
expression in the 2LYS case (Yq + Y() 2 , it is really only justified to compare the dominant 
terms Y Y and (YJ) 2 . 

In Table [TX] we show the corresponding fits to the matter distribution giving mo,2LYS = 
0.63 ± 0.05 and wii^lys = 1-0 ± 0.2 GeV. // these are now interpreted as the masses of 
constituent quarks, then scalar mesons constructed from two such non-interacting quarks 

- as in the naive quark model — would predict = 1.26 ± 0.10 and m{ = 2.0 ± 0.4 GeV. 
The latter is consistent with m^(2LY) = 2.0 ±0.3 GeV discussed earlier. However, the value 
of 777.Q is distinctly smaller than the earlier 7?7,q(2IY) = 1.53 ± 0.09 GeV and supports the 
point of view that scalar mesons cannot be described by the naive constituent quark model. 

4) GS and LGS. The above Yukawa and exponential forms arise naturally in quantum 
field theory, whereas gaussians do not. However, when — as in Ref. [[| - an attempt is 
made to understand the densities in terms of solutions of the Dirac equation, it would be 
difficult to reconcile an exponential or Yukawa asymptotic tail with the usual form of linearly 
rising confining potential cr. When such a potential is introduced as a scalar potential, the 
solutions of the Dirac equation are asymptotically gaussian. This is most easily seen when 

- in the notation of Ref. JTH] — the coupled Dirac equation for large r is written as 



- m(r)G(r) = -F'(r) and m(r)F(r) = G'(r) giving G" + (cr) 2 G = 0, (20) 

where m(r) = m + cr —>■ cr . The functions G and F are, therefore, seen to decay asymptot- 
ically as gaussians. Of course, the concept of a linearly rising confining potential does not 
hold for sufficiently large r, since eventually this will be quenched by the creation of qq pairs. 
Unfortunately, the actual demonstration of this unavoidable effect has yet to be achieved in 
a completely convincing manner for full QCD. However, the indications are that this will 
only occur for some value of r greater than about 1.2 fm, which — with the present lattice 
spacing of a ~ 0.14 fm — corresponds to a distance of almost 10a. Therefore, in the range 
of interest here (4a < r < 7a) the linearly rising potential is still expected to be important 
and its repulsion could well suppress the density from being an exponential decay to more 
like a gaussian decay. 

Support for a Dirac equation description is also given by our result that the charge and 
matter distributions are different — a feature not easy to understand in a non-relativistic 
approach. In the notation of Eq. ^ the charge and matter distributions can be expressed 
as x^{r) = G a {r)Gp{r) + F a (r)Fp(r) and x£f ( r ) = G a (r)Gp(r) — F a (r)Fp(r) respectively. 
Attempts are now underway to study to what extent the above distributions can indeed be 
interpreted in terms of solutions (G, F) of the Dirac equation. 

Our use, in section IT I [TJ, of the separable approximation when analysing the data was 



suggested by the form of the three-point correlation function C3 and supported by the 
presence in the data of the symmetry R(T) = R(T + 1) when T is odd — see Eq. [TO . 



However, as seen in Table [TJ this support is mainly for small values of r and is not present for 
the larger values of r. The above interpretation of the data in terms of Dirac wavefunctions 
is consistent with this — namely — for small r the lattice data gives x 1 ^ ~ 2^. This 
implies that G\{r) ^> F x (r) and both densities are approximately described by Gi{r)Gi{r) 

- a separable form. Whereas at the largest values of r studied here, we get x]} ~ 3x^ - 
implying that G\{r) ~ 2F 1 2 (r) and so destroying the simple separability of the densities. 
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D. More on the charge sum rule. 



Comparing the values of the sum rule in Table [V], where the estimates are made by 
directly summing over all the lattice as in Section [III Cj , and the values of J F in Table |VT|, 
where they are made by summing the separate contributions from each vertex (x, y, z), 
we see that the latter in the 2LY case can be as large as 1.7(3) — a number that seems to 
be slightly larger than the direct sum of about 1.4(1). Of course, within the quoted error 
bars these two estimates agree. Even so there are reasons — to be discussed below — why 
exact agreement is not necessary. Firstly, it must be remembered that the direct estimate 
includes contributions from the whole lattice i.e. upto values of r = 8\/3, whereas the data 
for a given r is based on only 18 values extending upto almost r = 6. More explicitly, the 
fits using Eqs. [T4] and |T5] are to data with r < 6 — the largest being x = 5, y = 3, z = 0. 
In fact, a more correct statement would be that the fits are mainly dictated by the data for 
r < 4, with the remaining data, which have relatively large error bars, playing more of a 
supportive than a decisive role. Therefore, there is no guarantee that the expressions based 
on Eqs. [14| and are a good representation beyond r pa 4. 

To test the importance of the density contributions from large values of r, the summation 
over r in Eq. [17| is carried out explicitly but truncated in two ways by introducing either a 
cubic cut-off (C) or a spherical cut-off (S). In the cubic case, the x, y and z sums are each 
limited to the range — r c to + r c , where r c takes on values ranging from 3 to 7. The first 
values of r c cover much of the range over which Eqs. |14| and [15] were used to fit the data, 
whereas by going to r c = 7 additional points are included that are outside this range. In the 
spherical cut-off, only values of x, y and z with x 2 + y 2 + z 2 < r 2 are kept in the summation. 
The outcome as seen in Table XI shows several points: 

1) There is little difference between the use of two lattice exponentials (2LE) versus two 
lattice Yukawas (2LY). This shows that, for the sum rule, the two forms are not only very 
similar at those values of (x, y, z) for which there exists lattice data but also at all other 
points on the lattice. Of course, when the expressions in Eqs. [b| and [15] are used over the 



whole lattice, there is no need to resort to the explicit summation in Eq. [18], since the result 
(Eq. [19]) is known. However, it did serve as a numerical check. 



2) The effect of including — where ever possible from Table [0 I - the 18 actual lattice 
data points instead of the fitted forms has the minor effect of decreasing the sums by about 
0.02. This is yet another reflection that the fitted forms are a good representation of the 
lattice data. Of course, these 18 points are only a small fraction of the total data needed 
in Eq. [17|. Fortunately, about one half of the sum rule comes from contributions within a 
volume where r < 3 and these have all been measured directly. 

3) The effect of using the Cubic (C)- versus the Spherical (S)-cutoff is large, with the 
latter being consistently about 0.2 smaller. This is natural, since for the same value of r c 
the C-cutoff embraces more lattice points. 

4) However, the most significant point is that the sums continue to increase significantly 
as r c goes beyond the range where the lattice data is measured — with the value at r c = 5 
for the S(C)-cutoff being about 1.1(1.4) increasing to 1.4(1.6) at r c = 7. 

The outcome is that, from the region r < 4, where the fits in Eqs. |14| and [15] are most 
reliable, the contribution to the sum rule is about 1.2 - 1.4, which already is consistent 
with the direct summation estimate in Table |V|. This, therefore, means that there is a 
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significant contribution of about 0.4 - 0.5 coming from values of r greater than 4. In detail, 
a contribution of about 0.2 comes from the range 4 < r < 5 and about 0.1 from the range 
5 < r < 6, where the fits seem to be supported by the data. The remaining discrepancy of 
almost 0.2 then comes from the periphery with r > 6, where there is no data to check the 
meaningfulness of the fits. 

It, therefore, seems that the direct sum rule in Table [V] could well be slightly smaller 
that the explicit sum in Eq. [T7|. This possible difference can be interpreted in two ways: 
The explicit sum is an overestimate or the direct sum is an underestimate. 

1) The fits using the expressions in Eqs. 14" and |T5| may indeed give a good fit to the data 
upto r « 4, but are not good estimates of the poor and missing data for r > 4. To test this, 
the lattice forms 2LE and 2LY are used upto r = r c and the single non-lattice forms based 
on 2E, 2Y and 2G in Table |VT] are used for r > r c . The outcome is shown in Table |X11 . 
As expected the combinations 2LE+E and 2LY+Y give results in the range 1.4-1.5, which 
are close to those of the spherical cutoff with r c = 7 in Table XI. On the other hand, the 
combinations 2LE+G and 2LY+G — containing the gaussian tail — give a result of ~ 1.3, 
which is noticeably smaller than the 2LE+E or 2LY+Y numbers. Furthermore, all of these 
combinations result in sumrules that are smaller than the r c = 7 cubic cutoff values of almost 
1.7. 

2) The second — and less likely - - interpretation is that the fits with Yukawa and 
exponential forms do indeed give a good estimate of the poor and missing data for r > 4. 
This would mean that the values of the sum rules in Table [V| are an underestimate and 
that the measured contributions there from r > 4 are too small. An example of this will be 
seen later in Table [Xllfj of the next subsection. There the three-point correlation function 
for r = 5 is consistent with that for (x = 3, y = 4) for T < 6 — but with errors that are 
twice as big. However, for T > 6 — the T range necessary for extracting densities — the 
signals from the two cases differ greatly, with those for r = 5 simply disappearing. Possibly 
similar underestimates could occur for larger values of r, where the fits from Eqs. |T4] and ITS] 
suggest a significant contribution 0.3) to the sumrule should arise. However, this would 
be surprising since it is not usual for there to be such systematic trends. Normally in such 
cases, one would expect the results to fluctuate from being too low for some some values of 
r to being too high at others. However, if the trend suggested by the r = 5 data were true 
then this would mean that the charge sumrule could be larger than that measured directly 

- possibly upto about 1.7. Unfortunately, to now get the value of unity expected in the 
continuum limit, would then require a vertex renormalisation factor of about 1/1.7 ~ 0.6 - 
a number significantly smaller than the estimates in Ref. |JTT . 



The first interpretation has two nice features. Firstly, with a charge sumrule of about 
1.3, the required vertex renormalisation factor to ensure unity in the continuum would be 



about 0.8 — a value more in line with the estimates in Ref. [[Ll|. Secondly, when — as 
in Ref. and discussed above — an attempt is made to understand the charge density 
in terms of solutions of the Dirac equation, exponential and Yukawa forms do not arise 
naturally, whereas gaussian forms do as illustrated in Eq. 

For the matter sum rule we saw in Table [V| that 0.9(1) was a reasonable compromise. 
The predictions using the algebraic fits are shown as If in Table |V11| and are seen to have 
much larger error bars than their charge counterparts. All that can be said is that the J F 's 
are — within these large uncertainties — consistent with the direct measurement of 0.9(1). 
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The conclusion from this section is not definite. If the asymptotic form of the density 
is indeed exponential or Yukawa, then the charge sumrule is 1.4 to 1.5. However, it is not 
possible to rule out an asymptotic form that is gaussian. In that case the charge sumrule 
could be less than 1.3. 



E. Rotational invariance 

Since off-axis points are considered, in principle it should be possible to check rotational 
invariance. In particular, comparisons between the on-axis point r = 5 and the off-axis point 
(x,y) = (3,4) data and also the r = 3 and (x,y,z) = (2,2, 1) data are of special interest. 
In Table |XIII| the ratios (C 3t p 1 F 1 (T, r) )/ (C*2,FiFi (T)) are shown for the dominant charge 
correlation. This table shows several points: 

1) The r = 3 and (2, 2, 1) data both have good signals for all values of T. This 
shows clearly that rotational invariance is violated with x^ iFi (r = 3) being almost twice 
x FiFi(2,2, 1). The reason for this lack of rotational invariance in such an extreme off-axis 
case could be due to the presence of the single z = 1 step in (2, 2, 1). A similar effect can 
be seen in Figure |3| for the (1, 1, 1) data which is also lower than the general trend. For less 
extreme cases such as (2, 1), rotational invariance seems to be better satisfied. 

2) The r = 5 signal disappears at T = 7 - rendering it useless to extract the density, since 
this requires reliable data for T > 7. This negative result is anyhow of interest, because in 
Ref. r = 5 was the largest value of r analysed and it suggested that the charge density 
was considerably smaller than would be expected from a simple exponential dependence. 
This is no longer the case. 

3) The (x, y) = (3, 4) signal is good for all T. This improvement over the r = 5 case 
is presumably due to the fact that each off-axis correlation is measured 24 times compared 
with the 6 on-axis measurements. 

4) For T < 7, the two sets of data are in agreement, as would be expected if rotational 
invariance had been achieved. In fact, this comparison can be pushed further by noting, 
that at these 3 values of T, C 3 (r = 5) appears to be slightly larger than C 3t -F 1 -p 1 (3,4:) — as 
would be suggested by the strong coupling model. 



If the lattice Yukawa and exponential forms in section IV are used to estimate the 



rotational invariance, we get the results in Table [XI V|. Both forms give the same results 



well within 1%. However, for the two cases considered involving directly the lattice data - 
r = 3 versus (2, 2, 1) and r = 5 versus (3, 4) — the off-axis values are about 10% less than 
the corresponding on-axis value. For r = 3 versus (2, 2, 1) this difference is considerably 
smaller than the actual lattice data requires. On the other hand, for r = 5 versus (3, 4) the 
reverse could be true with the difference being much larger than is suggested by the small 



T values in Table XIII 



It should be added that this study of the (2, 2, 1) case was only carried out after the 
bulk of this work was completed. In particular, the (2,2,1) data was not included in the 
fitting that led to the exponential, Yukawa and Gaussian parameters in Tables [VT| to [TX . 



However, its direct inclusion did not result in values of the a; and that were outside of 
the error bars given in Tables [VT] to |TX], eventhough each x 2 / n dot did increase by about 0.3 



- 0.4. A more realistic way to treat this data is to combine the r = 3 data with that for 
(2,2,1) using the relative weights of 6 to 24 — the number of equivalent possibilities for the 
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two cases. When this weighting is also used for the fitting functions, the x 2 / n dof decreases 
slightly from the values in Tables |VI] to fX[ However, again the final values of the a, and fj 
are unchanged within the errors quoted. This suggests that the overall fits are sufficiently 
stable that the effect of any particular case may lead to a large x 2 contribution for that case 
but still leave the fit essentially unchanged. It should be noted that the r = 3 and (2, 2, 1) 
charge data only differ by about two standard deviations, so that the above problem could 
even evaporate if more gauge configurations were used. 



V. CONCLUSION 

In this paper charge (vector) and matter (scalar) radial distributions have been measured 
on a lattice for the heavy-light meson (Qq), where Q is a static quark and q has a mass 
approximately that of the strange quark. The charge distribution could be determined 
reasonably well upto an interquark distance of about 6 lattice spacings i.e. ~ 0.8 fm. In 
comparison, the matter distribution measurements could only be carried out upto about 4 
lattice spacings i.e. « 0.6 fm. The drop-off of the charge distribution can be well described 
by the exchange of a vector meson of mass m 1 GeV. On the other hand, the drop-off of the 
matter distribution is described by the exchange of an scalar meson of mass « 1.5 GeV. 

In the conclusion of Ref. several refinements and extensions to that pilot calculation 
of charge and matter distributions were listed. Here we have carried out a few of these: 

1) Probably the most important "refinement" is the replacement of the quenched ap- 
proximation by the use of dynamical fermions. However, as seen in Figure ^|, we find that 
the two appear to be indistinguishable within the accuracy of the present work. In Refs. M 



and |T(J it was suggested that in the matter distribution there could be a difference due to 
the presence of disconnected quark-loop contributions. The fact that this is not seen here 
could be due to our use of sea quarks that have a mass about that of the strange quark 
mass. 

2) Radial correlations at off-axis points are now measured. This meant that the number 
of data points accessible before the noise takes over is much larger — going from about 6 to 
18. This enables us to achieve better algebraic fits to the data. However, one of our hopes, 
to see rotational invariance by comparing the data for r = 5 and (3,4) was only partially 
successful — see Table [Xlll| . 

3) The lattice spacing is now smaller — 0.14 fm compared with the earlier 0.17 fm. Also 
the number of gauge configurations is larger — 78 versus 20 earlier. 

However, the list in Ref. || contained other points not touched here: 
1) So far we have only extracted S'-wave correlations. We still need to measure the 
P1/2, P3/2, D 3 / 2 , -D5/2, ... densities corresponding to the energies extracted in Ref. ||. Also 
for a given orbital angular momentum, do these correlations show the degeneracy predicted 
in Ref. HI? 



2) The measurement of correlations in the baryonic and (Q 2 q 2 ) systems. Are these 
similar to those in the (Qq) system - as is the case when comparing correlations in few- 
nucleon systems? If this is so, then it would encourage phenomenological approaches such 
as that mentioned above utilizing the Dirac equation. 

3) As with all lattice calculations, there is the need to check the continuum limit by using 
finer and larger lattices. The former has, to some extent, been checked by the comparison 
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between our earlier work in Ref. with a = 0.17 fm and this study with a = 0.14 fm. 
There, as seen in Figure ||, the two sets of results agree after scaling. However, the two 
calculations also differ by their use of quenched versus unquenched gauge configurations. 
A more correct comparison would involve the use of the same type of configurations - 
quenched or unquenched. Even so, our present comparison — in spite of its failings — is 
encouraging. 

4) So far we have dealt with light quarks (valence and sea) that have the isospin form 
of u, (i-quarks but with a mass about that of the strange quark and heavy quarks that are 
static. The use of such quark masses means that the nearest physical meson with which we 
can possibly compare our calculations is the _B S (5.37 GeV). However, for a more realistic 
comparison we should eventually have a heavy quark with mass about 5 GeV and u, (i-quarks 
with their correct mass. 

5) In this work we have only probed the charge and matter distributions using the 74 and 
unity probes. However, other probes are possible such as: i) The pseudo-vector operator 
(7^75) needed for the B*Bn coupling — see Ref. |18| ; ii) Probes to study the color structure 



of the gluon fields and possible qq condensates surrounding the Qq system. 

This study is now at the stage where we have the energies and the corresponding charge 
and matter radial distributions for the ground and first excited S-wave states. In the future 
we hope to have these quantities for the other partial waves in Ref. [Q. These results can 
possibly be utilized in at least two ways: 

1) When, for example, calculating electromagnetic transitions between different Qq states 
the form of the transition matrix elements could be guided by the above radial distributions. 

2) As discussed in the Introduction — and also as partial motivation for the separable 
method for analysing the lattice data — the radial distributions could be interpreted in terms 
of wavef unctions. This would mean we are in the position of having some of the eigenvalues 
and eigenf unctions of an Hamiltonian, whose form we do not know. We are then at liberty 
to find this Hamiltonian to construct an effective theory. This effective Hamiltonian could 
be, for example, of Schroedinger or Dirac form with suitable interactions. As discussed in 
Section [IV Q , the latter form is probably more appropriate, since the charge and matter 
distributions are different — a feature not easy to understand in a non-relativistic approach. 

This study has shown that — using lattice techniques — reliable estimates can be made 
not only of spectra but also of wavefunction information. This is just the beginning - 
with future studies expected to enlarge this information and also to attempt interpretations 
outside of quantum field theory. 
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APPENDIX A: THE SYMMETRY R(T) = R(T + 1) FOR ODD T 

In the separable case where x al3 (r) = y a { r )yp{ r )i for even T 
R(T, r) = - In 

+ ln 1 + (r)° ' ■ * T c(l,l) yi (ry 



I I c(2,l)y 2 (r) ( m2 - mi )^ . c(3,l)y 3 (r) -( m3 - mi )% 

L c(l,l)yi(r) ~ c(l,l)j/i(r) 

J + c(2,l)y 2 (r) c -( m ^- mi )^ + e (3,l)y 3 ( r ) c _( m3 _ mi )2^2 



whereas for odd T 
R(T,r) = 



-ln 
+ ln 



I + e(2,l)i/a(r) c _( m2 - TOl + £ftlkM e -K-mi)^ 



c(l,l)yi(r) 



c(l,l)j/i(0 



J + c(2,% 2 (r) c -( m ,- mi )3^1 + c(3,l)y 3 (r) p -(m 3 -mi )^ 



c(l,l)yi(r) 



c(l,l)3/iW 



(Al) 

(A2) 

(A3) 
(A4) 



Here we have dropped the fuzzing indices on R. By inspection it is seen R(T) = R(T + 1) 
as given in Eq. [lOj and demonstrated in Table ||. 
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TABLES 



TABLE I. Values of the parameters am a and vf, needed for fitting the (^-correlations. 
Case Q: 20 quenched configurations with Mi = 3 and T2 jmm = 4 — referred to as Case 3 in || 
Case A: 78 dynamical configurations with Mi = 3 and T2 mi n = 4 
Case B: 78 dynamical configurations with Mi = 4 and T2 mi n = 3 
Case C: 78 dynamical configurations with Mi = 4 and Ti :Ui \ n = 4 

The entries {...} show the values when correlations are removed as discussed in the text. 
The entries [..] for Am2i and A77131 include the ratio of lattice spacings 0.14/0.17. 
The entries marked with a dash are not applicable for Mi = 3. 



am a 


Case Q 


Case A 


Case B 


Case C 


vf 


3x3 T m i n = 4 


3x3 T m ; n = 4 


4 x 4 T min = 3 


4 x 4 T min = 4 


ami 


0.8721(19) 


0.8580(11) 


0.8340(40) 


0.8280(88) 


ami 


1.263(13) 


1.2267(51) 


1.166(11) 


1.138(35) 


am% 


1.94(30) 


1.93(13) 


1.632(42) 


1.52(12) 


am^ 
vi 

Fi 


0.4847(56) 
1.519(10) 
0.8402(38) 


0.4344(31) 
1.3779(70) 
0.8008(21) 


1.889(58) 
0.3757(99) 
1.227(26) 
0.731(12) 


1.85(18) 
0.359(24) 
1.181(65) 
0.711(31) 


L 


0.816(16) 


0.8405(65) 


0.801(12) 


0.757(55) 


V 2 

vl 

4t 


0.644(49) 
-0.251(33) 
-0.28(22) 
2.2(1.4) 


0.874(22) 
-0.115(12) 
-0.348(97) 

2.34(63) 


1.185(44) 
0.169(47) 
-0.459(29) 
1.36(16) 


1.260(83) 
0.26(12) 
-0.534(85) 
0.86(53) 


4 2 

vf 


-1.13(81) 


-0.84(39) 


0.56(19) 
0.000(69) 
-0.45(33) 

1.89(22) 


0.61(25) 
0.05(19) 
-0.73(70) 
1.64(78) 


"2,data 


54 


48 


54 


48 


^2, par am 
n 2,do{ 


12 
42 


12 

36 


16{7} 
38{45} 


16 
32 


X 2 /n 2 ,doi 


0.65 


4.27 


0.16{0.13} 


0.14 


Amu 


0.391(13) 
[0.322(H)] 


0.369(5) 


0.332(12) 


0.310(36) 


Am 3 i 
Am4i 


1.07(30) 
[0.88(24)] 


1.07(13) 


0.798(42) 
1.05(6) 


0.69(12) 
1.02(18) 
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TABLE II. The ratio R(T) defined by Eq. | and having the symmetry R(T) = R(T + 1) in 
Eq. |l0], when T is odd. Here R(T) is calculated for C^FiFi — the dominant correlation. 



T 


R(T,r = 1) 


R(T,x = l,y = l) 


R(T,x = l,y = l,z = l) 


R(T,x = 3,y = A) 


5 


1.0222 


0.958 


0.890 


0.611 


6 


1.0223 


0.955 


0.882 


0.585 


7 


0.9536 


0.919 


0.876 


0.598 


8 


0.9507 


0.913 


0.874 


0.544 


9 


0.9166 


0.879 


0.844 


0.590 


10 


0.9133 


0.857 


0.839 


0.524 



TABLE III. Estimates of the charge radial distributions [x 11 (r)] for the ground state. 
Separable case (S): M 2 = 4 in C 2 , M 3 = 3 in C(3). T 2imin = 3 and T 3jmin = 6, 8. 
Separable case (S): M 2 = 4 in C 2 , M3 = 3 in C(3). T 2im i n = 4 and T3 jm ; n = 8. 
Non-separable case (NS): M 2 = M3 = 4 in C 2 and C3 but x 13 = 0, when neither i nor j is 1. 
22 im i n = 3 and Ts jm - m = 6, 8. 

The column labelled Best Estimate is a summary of the previous 5 columns. 

The column labelled Approx. Estimate is from Eq. [l2| for M 2 = 3 — Case A in Table |[ 

The line marked as 'ave' is the weighted mean [6x n (3, 0, 0) + 24x 11 (2, 2, l)]/30 











s 


S 


S 


NS 


NS 


Best 


Approx. 


^2,min 




3 


3 


4 


3 


3 


Estimate 


Estimate 


^3,min 




6 


8 


8 


6 


8 






r 


X 


y 


z 








0.00 











0.0378(2) 


0.0377(4) 


0.0349(5) 


0.0428(3) 


0.0401(8) 


0.038(4) 


0.047(2) 


1.00 


1 








0.01214(10) 


0.0127(2) 


0.0112(3) 


0.01343(15) 


0.0135(4) 


0.0125(15) 


0.015(1) 


1.41 


1 


1 





0.00704(8) 


0.0075(2) 


0.0068(2) 


0.00831(10) 


0.0091(3) 


0.0080(15) 


0.009(2) 


1.73 


1 


1 


1 


0.00501(8) 


0.0052(2) 


0.0048(2) 


0.00594(11) 


0.0063(3) 


0.0055(7) 


0.0060(3) 


2.00 


2 








0.00749(9) 


0.0072(2) 


0.0070(3) 


0.00627(11) 


0.0068(3) 


0.0070(5) 


0.0071(5) 


2.24 


2 


1 





0.00476(6) 


0.00470(13) 


0.0045(2) 


0.00464(6) 


0.0047(2) 


0.0047(2) 


0.0046(3) 


2.83 


2 


2 





0.00307(7) 


0.0031(2) 


0.0031(2) 


0.00299(7) 


0.0031(2) 


0.0032(2) 


0.0028(2) 


3.00 


3 








0.00329(10) 


0.0034(3) 


0.0035(3) 


0.00295(11) 


0.0035(3) 


0.0034(4) 


0.0032(5) 


3.00 


2 


2 


1 


0.00221(6) 


0.00216(15) 


0.0021(2) 


0.00215(6) 


0.0019(2) 


0.00215(15) 


0.0021(2) 


3.00 


ave 




0.00243(6) 


0.00241(13) 


0.0023(2) 


0.00231(5) 


0.0022(2) 


0.00240(14) 


0.0023(2) 


3.16 


3 


1 





0.00234(6) 


0.00244(14) 


0.0025(2) 


0.00226(6) 


0.0026(2) 


0.0025(3) 


0.0024(2) 


3.61 


3 


2 





0.00164(6) 


0.00178(13) 


0.0019(2) 


0.00166(5) 


0.00188(14) 


0.0018(2) 


0.0017(2) 


4.00 


4 








0.00141(10) 


0.0015(3) 


0.0016(4) 


0.00126(10) 


0.0014(3) 


0.0015(3) 


0.0012(2) 


4.12 


4 


1 





0.00095(5) 


0.00105(13) 


0.0011(2) 


0.00092(5) 


0.00115(15) 


0.0011(2) 


0.0009(1) 


4.24 


3 


3 





0.00082(8) 


0.0010(2) 


0.0010(3) 


0.00085(7) 


0.0011(2) 


0.0010(3) 


0.0008(2) 


4.47 


4 


2 





0.00074(5) 


0.00075(12) 


0.0007(2) 


0.00076(5) 


0.00068(15) 


0.00070(15) 


0.0005(2) 


5.00 


4 


3 





0.00067(6) 


0.00080(10) 


0.0009(2) 


0.00078(5) 


0.00115(14) 


0.0009(3) 


0.0008(2) 


5.10 


5 


1 





0.00046(6) 


0.00053(14) 


0.0006(2) 


0.00050(5) 


0.0007(2) 


0.0006(2) 


0.0005(2) 


5.39 


5 


2 





0.00033(8) 


0.00044(14) 


0.0005(2) 


0.00036(5) 


0.0006(2) 


0.0005(2) 


0.0003(1) 


5.83 


5 


3 





0.00020(5) 


0.00026(12) 


0.0003(2) 


0.00027(5) 


0.0005(2) 


0.00035(15) 


0.0004(2) 
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TABLE IV. Best estimates of radial distributions of charge and matter distributions involving 
also excited states. The column labelled x 11 (AE) is an approximate estimate from Eq. [TJ for 
M2 = 3 — Case A in Table |[ The entry '-' implies that a reasonable signal could not be obtained 
and the entry 'ave' has the same definition as in Table n. 



Charge 


Matter 


r 


X 12 


x 13 


X 22 




x 11 (AE) 


x 12 


0.00 


0.079(4) 


-0.050(2) 


0.18(2) 


0.045(3) 


0.058(2) 


0.090(5) 


1.00 


0.0160(5) 


0.015(2) 


0.021(3) 


0.0133(5) 


0.016(1) 


0.009(1) 


1.41 


0.0045(15) 


0.010(2) 


0.003(1) 


0.0073(7) 


0.0076(6) 


-0.001(1) 


1.73 


0.0010(10) 


0.007(2) 


0.005(4) 


0.0049(5) 


0.0045(6) 


-0.0030(7) 


2.00 


0.0020(5) 


0.0005(15) 


0.0005(4) 


0.006(1) 


0.0054(5) 


-0.003(1) 


2.24 


0.0000(3) 


0.0022(5) 


0.0000(1) 


0.0038(4) 


0.0032(3) 


-0.0030(7) 


2.83 


-0.0010(3) 


-0.0005(10) 


0.0003(2) 


0.0017(2) 


0.0013(3) 


-0.0015(3) 


3.00 


-0.0009(4) 


-0.0018(15) 


0.0003(2) 


0.0014(3) 


0.0013(3) 


-0.0010(3) 


3.00 


-0.0007(2) 


0.0010(10) 


0.0003(2) 


0.0016(2) 


0.0009(2) 


-0.0018(5) 


ave 


-0.0007(2) 


0.0004(9) 


0.0003(2) 


0.0016(2) 


0.0010(2) 


-0.0016(4) 


3.16 


-0.0009(2) 


-0.0010(10) 


0.0003(2) 


0.0009(2) 


0.0009(3) 


-0.0005(5) 


3.61 


-0.0010(2) 


-0.0010(5) 


0.0006(2) 


0.0007(2) 


0.0005(3) 


-0.0008(4) 


4.00 


-0.0007(2) 


-0.0018(12) 


0.0005(3) 








4.12 


-0.00045(15) 


-0.0010(8) 


0.0003(2) 








4.24 


-0.00050(15) 


-0.0011(8) 


0.0005(3) 








4.47 


-0.00040(15) 


-0.0000(5) 


0.00025(15) 








5.00 


-0.0007(2) 


-0.0010(5) 


0.00065(15) 








5.10 


-0.00040(15) 


-0.0008(4) 


0.0004(2) 








5.39 


-0.0003(2) 


-0.0005(3) 


0.0003(2) 








5.83 


-0.0003(2) 


0.0001(3) 


0.0001(2) 
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TABLE V. The charge and matter sum rules in the ground state x 11 and between excited 
states. 



The symbols N and NS are as in Table E3 . 

Q refers to the quenched results (case 3) in Ref. ||. 

The entries marked with a dash are not applicable for M3 = 3. 



T-1 , min / -^3 , min 

x al3 


S S S 
(3/6) (3/8) (4/8) 


NS NS NS Q-NS 
(3/6) (3/7) (3/8) (3/8) 


Charge 






x 12 
x 13 
x 14 
x 22 


1.12(3) 1.26(6) 1.26(8) 
0.09(3) -0.10(9) -0.05(12) 
0.02(9) -0.3(4) -0.1(4) 

0.008(5) 0.007(14) 0.002(9) 


1.29(4) 1.36(6) 1.42(11) 1.41(5) 
-0.28(7) -0.40(12) -0.5(2) -0.5(2) 
0.09(11) -0.2(2) -0.3(6) 0(1) 
-0 1(2) —0 4(4) -0 5(1 1) 
0.8(2) 0.9(3) 1.1(6) 0.9(9) 


X 2 / n dof 


2.52 0.65 0.63 


0.31 0.13 0.09 0.26 


Matter 






X 11 
X 12 
X 13 

x 14 

X 22 


0.66(5) 0.85(14) 1.0(2) 
-0.15(6) -0.5(2) -0.6(3) 
-0.4(2) -1.7(1.0) -1.6(9) 

0.04(3) 0.3(2) 0.4(4) 


0.73(8) 0.83(11) 1.1(2) 0.38(15) 
-0.27(13) -0.4(2) -0.8(4) -0.1(5) 
-0.3(2) -0.9(4) -2.7(1.2) -1(5) 
-0.2(3) -0.8(7) -3.2(2.2) - 
0.2(3) 0.2(6) 0.2(1.2) 1(2) 


X 2 /ndoi 


0.46 0.29 0.25 


0.31 0.25 0.12 0.35 



TABLE VI. Fits to the charge lattice data with exponential, Yukawa and Gaussian forms. 
Column 2E refer to the exponential form ^ 1 af exp(— r/rf ). 
Column 2Y refer to the Yukawa form ^ 1 aj ex.p(—r/rj)/r. 
Column 2G refer to the gaussian form ^ 1 af exp(— (r/r^) 2 ). 

Also fits to the lattice data with the lattice Yukawa (2LY), exponential (2LE) Gaussian (2LG) 
forms from Eqs. [l4|, ^ and [l6|. 

J E = 8vrEo,i af{rf)\ I Y = 4^Eo,i a T( r 7) 2 ' *S = ^ Eo,l «?K°) 3 



are the spacial integrals of these functions. 

The entries marked fixed in the minimization. 



Form(i ? ) 


2E 


2LE 


2Y 


2LY 


2G 


2LG 


a l 


0.0250(22) 


0.0245(21) 


0.069(10) 


0.066(3) 


0.0086(5) 


0.0075(7) 


rl 


1.36(6) 


1.37(6) 


1.70(6) 


1.78(7) 


2.92(9) 


3.07(12) 


a\ 


0.013(5) 


0.373* 


-0.073* 


-0.060* 


0.0293(41) 


0.019(4) 


F 
r\ 


0.2(4) 


0.116(15) 


1.00(17) 


1.12(11) 


0.75(7) 


0.99(15) 


I F 


1.6(3) 


1.6(3) 


1.6(5) 


1.7(3) 


1.3(1) 


1.3(2) 


X 2 / n dof 


1.32 


0.81 


1.33 


0.94 


1.47 


0.93 
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TABLE VII. Fits to the matter lattice data with exponential, Yukawa and Gaussian forms. In 
2LY^ the r = data point is not fitted. Other notation as in Table VI. 



Form (F) 


2E 


2LE 


2Y 


2LY 1 " 


2G 


2LG 


«o 


0.036(20) 


0.0345(26) 


0.178(12) 


0.186(23) 


0.0105(12) 


0.0101(16) 


F 
r 


0.91(13) 


0.938(39) 


1.04(12) 


0.92(6) 


2.11(9) 


2.14(11) 


a F 


0.009(20) 


0.561* 


-0.168* 


-0.207* 


0.0346(32) 


0.034(10) 


rf 


0.45(55) 


0.099(11) 


0.90(18) 


0.71(12) 


0.71(4) 


0.70(12) 


I F 


0.7(5) 


0.73(10) 


0.7(9) 


0.7(6) 


0.62(10) 


0.62(13) 


X 2 / n dof 


1.80 


1.16 


1.74 


1.11 


1.11 


1.00 



TABLE VIII. Fits to the separable function y\(r) defined as yi(r) = (x 11 ) 1 ^ 2 , where the x 
are the charge density values labelled as Best Estimate in p| . Notation as in Table [V|. 



Form (LFS) 


2ES 


2LES 


2YS 


2LYS 


2GS 


2LGS 


LFS 
"0 


0.159(7) 


0.156(8) 


0.47(3) 


2.4310(60) 


0.092(3) 


0.083(4) 


^LFS 
r 


2.73(12) 


2.59(9) 


4.40(21) 


2.86(11) 


4.19(12) 


4.44(17) 


a\ FS 


0.036(13) 


0.0422* 


-0.49* 


-2.4196* 


0.102(11) 


0.065(6) 


^LFS 


0.2(5) 


0.41(17) 


1.66(22) 


2.43(10) 


0.85(9) 


1.24(15) 


xV^dof 


1.32 


1.11 


1.33 


1.11 


1.75 


1.07 



TABLE IX. Fits to the separable function y\(r) defined as y\{r) 
the matter density values in IV. In 2LY^ the r 
Table (VT. 



(x 11 ) 1 / 2 , where the x 11 are 
data point is not fitted. Other notation as in 



Form (LFS) 


2YS 


2LYS f 


n LFS 


0.75(4) 


0.74(4) 


r LFS 


2.39(27) 


2.25(18) 


a\ FS 


-0.728* 


-0.774* 


^LFS 


1.52(28) 


1.36(20) 


X 2 /n<ioi 


1.84 


1.25 
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TABLE X. The value of rg from fits to the Fourier transform of the charge density x (r) as 



defined in Eq. 17 



Row A: ro directly from Tables VI and VII] . 

Row B: ro extracted from the Fourier Transform of x n (r) when these are expressed in terms of 
the analytic expressions in Eqs. [14|, |l5] and ^ with the parameters in Tables VI and |VTII . Option 
(i) is used for defining the x 2 . 
Row C: Same as Row A but using Option (ii). 

Row D: The x n (r) are the same as in row B and C, except that the data from Table III is used 

wherever possible — Option (i) used. 

Row E: Same as Row D but using Option (ii). 

The entry '-' means a sensible solution could not be found 



Extract with 




2LY 




2LE 


A 


1.78(7) 


1.37(6) 


Form 


2LY 




2LE 


2LE 


B 


1.78(24) 




1.75(24) 


1.38(22) 


C 


1.78(9) 




2.5(6) 


1.38(15) 


D 


1.78(20) 




1.75(22) 


1.39(23) 


E 


3.3(1.0) 






1.41(16) 



TABLE XL Estimates of the charge sum rule from the k = component of Eq. 17. The r 
summation is truncated in two ways — cubic(C) cut-off, where \x\, \y\ and \z\ are all < r c and 
the spherical(S) cut-off, where only values of x, y and z with x 2 + y 2 + z 2 < r 2 are kept in the 
summation. The columns labelled I 2E and I 2Y are the integrals from to r c of the functions 2E 
and 2Y in Table fVl. The row labelled None means the r summation is over the whole L 3 lattice. 







2LE+data 


2E 


2LY+data 


2Y 


r c 


C 


S 


j2E 


C 


S 


j2Y 


3 


0.95 


0.62 


0.63 


0.95 


0.62 


0.58 


4 


1.20 


0.87 


0.90 


1.20 


0.87 


0.85 


5 


1.37 


1.11 


1.10 


1.38 


1.11 


1.06 


6 


1.49 


1.30 


1.23 


1.51 


1.30 


1.21 


7 


1.58 


1.42 


1.31 


1.61 


1.42 


1.37 


None 


1.62 


1.62 


1.6(3) 


1.65 


1.65 


1.6(5) 
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TABLE XII. This is a continuation of Table [X| for the spherical cutoff case (S) containing the 
available data. Here single exponential (E), Yukawa (Y) and gaussians (G) tails are included for 
r > r c . The parameters of these tails are obtained from data with r > 3. The combined functions 
(2LE+E, 2LY+Y, 2LE+G and 2LY+G) are integrated using the cubic cutoff with r = 7. 



r c 


2LE+E 


2LY+Y 


2LE+G 


2LY+G 


3 


1.42 


1.46 


1.30 


1.30 


4 


1.40 


1.45 


1.28 


1.28 


5 


1.41 


1.46 


1.29 


1.29 





TABLE XIII. Check on 


the rotation invariance 


of R(r) = (C 3>FlFl (T, 


r))/(C 2)FlFl (T)) for 


r = 


= 3 


versus x = 2, y = 2, z 


= 1 and r = 5 versus x 


= 3, y = 4. All entries should be multiplied 


by 10 


-4 








T 




R(r = 3) 


R(x = 2,y = 2,z = l) 


R(r = 5) 


R(x = 3,y = A) 


4 




23.0(4) 


11.6(2) 


1.2(4) 


1.2(2) 


5 




24.4(4) 


13.4(3) 


1.7(4) 


1.6(2) 


6 




25.4(6) 


15.3(4) 


2.0(5) 


2.2(3) 


7 




25.7(7) 


16.2(4) 


1.1(6) 


2.8(3) 


8 




26.2(1.1) 


17.7(0.7) 


0.8(9) 


3.9(5) 


9 




26.8(1.5) 


17.4(0.9) 


1.0(1.2) 


5.0(7) 


10 




28.8(2.7) 


17.4(1.5) 


2.9(2.2) 


7.0(1.3) 



TABLE XIV. Lack of rotation invariance in the lattice Yukawa and exponential forms in Eqs. 14 
and 15 for r = 3 versus (x, y, z) = (2, 2, 1) and r = 5 versus (x, y) = (3, 4). 



Case 


r = 3 


(2,2,1) 


Ratio 


2LY(Table [V| 


0.002893 


0.002566 


0.89 


2LE (Table 0) 


0.002893 


0.002573 


0.89 




r = 5 


(3,4) 


Ratio 


2LY(Table [V| 


0.000709 


0.000653 


0.92 


2LE (Table 0) 


0.000706 


0.000651 


0.92 
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FIGURES 

FIG. 1. The four contributions to the two-point correlation function C2 

75 ** 75 75** 75 * 75 75* 




75 * 75 75 * 75 ** 75 75** 



a) b) c) d) 



FIG. 2. The two contributions to the three-point correlation function C3 
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FIG. 3. The ground state charge (C) and matter (M) densities [x n (r)] as a function of r/a. 
The lines shows a fit to these densities with a sum of two lattice exponential functions. The scaled 
quenched results of Ref. are also shown by filled circles and squares. 



Charge and matter densities x 




10 1 1 1 1 1 1 1 

1 2 3 4 5 6 

r 
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FIG. 4. a) and b) — the ratios x 12 /x u and x 13 /x u for the charge distribution, c) and d) these 
ratios for the matter distribution 
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FIG. 5. The separable " wavefunctions" y% and 1/2 denned in the density x a @{r) = y a (r)yf3(r). 

a) yi — Solid line for the charge, dashed for the matter, solid circles the Bethe-Salpeter wavefunc- 
tion from Ref. ||, 

b) y2 — notation as in a) 
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FIG. 6. Fit to the lattice data in Table III with lattice exponential (2LE), Yukawa (2LY) and 
gaussian (2LG) forms in Table VI. 




35 



